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by 

2)  2) 
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Dedicated  to  Eugene  Issacson  on  the  occasion  of  his  70th  birthday. 

The  problem  of  determining  closed  curves  that  are  invariant  under  zcrr.z 
given  mapping  arises  in  many  applications  and  various  numerical  techniques 
have  been  proposed  for  the  calculation  of  such  invariant  cycles.  We  refer  here 
only  to  Doedel  [1],  Iooss  et  al  [3],  Kevrekidis  et  al  {4],  Van  Veldhuizen ,  [ 6 ] , 
where  further  references  may  be  found. 

As  in  [4]  we  consider  a  continuously  differentiable  mapping 


(s,x)  e  S^xRm 


($(s)  ,  F(  s  ,  x)  )  ,  S> :  S 


FisW" 


(1) 


where  S1  is  the  unit  circle  in  Rm.  In  other  words,  (1)  is  a  mapping  of  the 
cylinder  S^xR™  into  itself.  Then 


7 :  S1  -*  S^xRm,  7<s)  -  (s,r(s)),  s  e  S^; 


(2) 


is  invariant  under  (1)  exactly  if 


r(#(s))  -  F(s , r(s) ) ,  s  e  S1  (3) 
As  usual,  for  the  computation  we  approximate  this  curve  by  some  polygon  with 
vertices  x.^  =  rCs^  ,  i  -  l,2,...,n.  But,  other  than  in  the  cited  references, 
we  also  discretize  the  circle  mapping  $  by  the  following  mapping  from 


1)  This  work  was  supported  in  part  by  the  Office  of  Naval  Research  under 
contract  N-00014-80-G-9455 ,  the  National  Science  Foundation  under  grant 
DCR- 8309926,  and  the  Air  Force  Office  of  Scientific  Research  under  grant 
84-0131. 

2)  Department  of  Mathematics  and  Statistics,  University  of  Pittsburgh, 
Pittsburgh,  PA  15260. 

......  .i/jt  let:/ 

Availability  Codes 
Avail  or, d/or 
jDlat  j  Special 


2 


For  the  numerical  solution  of  discretizations  of  the  problem  (3) , 

Newton's  method  is  frequently  employed.  But  then  the  computational  cost  is 

3 

typically  of  the  order  of  (mn)  and  hence  grows  quickly  with  the  number  of 
vertices  of  the  approximating  polygon.  Here  we  want  to  show  that  for  discreti¬ 
zations  of  the  form  (5)  it  is  possible  to  achieve  a  growth-rate  that  is  linear 
in  mn.  This  is  certainly  a  considerable  reduction,  especially  for  applications 
involving  periodic  solutions  of  large  systems  of  ordinary  differential  equa¬ 
tions  . 


Let  T  denote  the  connectivity  graph  representing  the  mapping  N  -  N  of 
(5);  that  is,  T  has  N  as  its  vertex  set  and  a  directed  arc  from  i  e  N  to 
j  e  N  exactly  if  j  -  t(i).  Since  t  is  a  mapping  of  N  into  itself,  the  out- 
degree  of  any  vertex  of  T  is  exactly  one,  and,  for  any  starting  point  i  e  N, 
the  directed  path  through  the  vertices 

'i-1'  k-°’  l----  (7) 
is  uniquely  determined.  Because  T  has  only  n  vertices,  any  such  path  ulti¬ 
mately  involves  repeated  vertices.  Suppose  that  the  first  q>0  vertices 
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i^,  k-0,...,q,  of  (7)  are  distinct  but  that  i^+^  repeats  one  of  these  earlier 
vertices.  Then  we  have  two  possible  cases,  namely, 

(i)  i  ,  -  i  ;  that  is  i  is  a  fixed  point  of  <t  which  implies  that  i,  -  i 

q+i  q  q  k  q 

for  k>q . 

(ii)  ”  ip  for  some  p,  0<p<q  .  Then  the  vertices  i^,  k-p , p+1 ,  .  .  .  ,  q  form  a 

non- degenerate ,  directed  cycle  of  T  . 

If  two  paths  (7)  with  distinct  starting  points  have  a  vertex  in  common, 
then  they  must  be  identical  from  then  on,  for  otherwise  there  would  have  to  be 
a  vertex  with  out-degree  two.  Hence,  in  particular,  no  two  distinct  cycles  can 
intersect.  Thus  the  vertices  i  of  T  belong  to  three  disjoint  classes: 


(a)  i  is  a  fixed  point  of  <t>  , 


(b)  i  belongs  to  a  unique  cycle  of  T  , 


(c)  there  is  a  unique  directed  path  through  i  which  terminates  either  at  a 
fixed  point  of  $  or  on  a  cycle  of  T. 

In  other  words,  T  consists  of  finitely  many  disjoint  components 
Tj ,  j-1 . r  each  of  which  contains  exactly  one  fixed  point  of  $  or  one  non¬ 

degenerate  cycle  of  T  and  all  other  vertices  of  I\  belong  to  the  class  (c) . 

If  Newton's  method  is  applied  to  the  system  (5),  then,  at  each  step,  a 
linear  system  of  the  form 


(T-A)  u  -  b  ,  u  -  (ux . un)T>  uieRm> 


Tu  -  («(Ul) . *(un))\ 


has  to  be  solved,  where 
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A  -  diag(A1 . An),  Al  -  DF(x.)  e  R  .  b  -(bj. . bn)i>  bt  -  F(x(i))  €  R® 

Evidently,  (8)  represents  a  recursion 

u  -  bi  +  Ai  ui  -  Lk+1  “  k“0'1’2 -  (9) 

k+1  k  k  k 

Thus,  for  any  vertex  i  of  T  of  class  (a);  that  is,  for  any  fixed  point  i  of  'l> , 

the  value  of  u.  can  be  calculated  from 
1 

(I  -  A. ;u.  -  b. 
m  i  i  l 

where  1^  is  the  identity  on  R®,  if  only  the  matrix  on  the  left  is  nonsingular. 

Let  i  be  some  vertex  of  T  of  class  (b)  or  (c).  Then  the  path  (7)  started  at 

ip  -  i  begins  with  q>0  distinct  vertices  i^,  k-0,l,...,q  and  we  have 

i  -  i  ,  0<p<q  . 
q+1  p 

We  define  recursively  the  quantities 


d  —  b.  +  A  d,  ,  B,  —  A.  B.  ,  k— 0 , 1 . q+1 

k+1  \  k  1k  ik+l  \  \ 

then  (9)  can  be  rewritten  in  terms  of  the  unknown  first  u-value  u^ 


u  -  d  +  B  u,  ,  k-0 , 1 , . . . q+1  (11) 

k  \  \  0 

Therefore,  u.  can  be  computed  from  the  equation 
t0 

(B.  -  B.  )  u,  -  d.  -  d.  (12) 

1  i  in  1  1 

p  q  0  q  p 

provided  that  the  difference  of  the  B-matrices  on  the  left  is  non-singular. 

Then,  in  turn,  all  u.  ,  k-l,2,...,q  can  be  calculated  from  (11).  Moreover,  for 

Lk 

any  other  vertex  in  the  same  component  of  T  as  ig  the  corresponding  path  will 
end  at  a  vertex  j  where  the  u-value  is  already  known.  Hence  from  (12)  we  can 


calculate  Uj  provided  that  Bj  is  invertible.  Thus,  by  continuing  this  process, 
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we  can  compute  the  u-values  at  all  vertices  of  T  as  long  as  no  singular 
matrices  are  encountered.  Since  this  approach  is  equivalent  with  a  recursive 
evaluation  of  the  solution  of  (8)  it  follows  readily  that  --in  real  arith¬ 
metic  --no  singular  matrix  will  ever  arise  if  the  matrix  T  -  A  in  (8)  is  non¬ 
singular  . 

Evidently,  this  recursive  solution  of  (8)  involves  only  the  solution  of 
linear  systems  of  the  order  m  and  hence  can  be  expected  to  be  considerably 
faster  than  any  standard  decomposition  of  the  full  matrix  in  (8).  In  fact,  all 
our  examples  show  that  the  run-time  of  the  process  is  approximately  linear  in 
run  . 


In  many  applications,  the  circle  map  $  in  (la/b)  is  a  fixed  rotation: 

s  -*■  s  +  a>  mod  2ir,  0  <  u>  <  2*  (13) 
Here  it  follows  directly  from  (4)  that  ♦  must  be  injective  and  hence 
represents  a  permutation  on  N  .  Thus  T  now  consists  of  a  finite  number  of  dis¬ 
joint  cycles  and  the  solution  of  (8)  involves  only  the  solution  of  finitely 
many  systems  of  the  form  (12) . 


As  a  first  sample  problem  for  this  case  we  consider  the  following  func¬ 
tion  (1)  used  also  in  [4] : 


1 

2  1 

$(s)  -  s  +  mod  2?r  ,  ca  -  2ff  (5  -1),  s  e  S  (14a) 

F(s,x)  -  X  x(l  -  x)  +  «  cos(s),  s  6  S1,  x  e  R1  (14b) 

where  X  -  3.46  and  e  is  a  parameter. 


A  code  was  implemented  on  a  VAX- 8650  which  incorporates  the  above  solu¬ 
tion  method.  It  uses  a  continuation  procedure  similar  to  that  described  in 


6 


[5]  to  obtain  invariant  cycles  for  a  sequence  of  c -values  starting  from  the 
constant  solution  x  -  0.44  for  «  -  0.001  .  Table  1  gives  some  of  the  computed 
values  of  r(0)  for  different  t  and  for  n  -  377 


r(0) 


r 

|  0.001 

l 

0.437523  | 

|  0.01 

0.411982  | 

|  0.02 

0.384902  | 

|  0.03 

0.376986  | 

|  0.032 

0.375482  | 

|  0.034 

0.373925  | 

|  0.035 

0.375158  | 

|  0.0351 

0.374768  | 

|  0.0352 

0.374059  | 

|  0.035245 

0.373697  | 

i 

TABLE  1 

All  calculation  used  double  precision  and  a  corrector  tolerance  of  10  ^  A 
corresponding  run  with  n  -  610  produced  values  of  r(0)  which  differed  only  by 
about  1.0x10  ^  .  The  elapsed  clock- time  T  for  the  execution  of  each  cycle  was 
recorded,  and  Table  2  shows  the  quotients  T/(nk)  for  different  dimensions  n 
where  k  denotes  the  number  of  Newton  corrector  steps  taken  by  the  process.  The 
linear  dependence  of  the  time  on  n  is  clearly  visible. 
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1  n 

233 

377 

610  | 

T 

T 

T  | 

1  € 

nk 

nk 

nk 

|  0.001 

2.58 

3.24 

4.36  | 

|  0.002 

2.58 

3.24 

4.36  | 

|  0.003 

2.58 

3.27 

4.37  | 

TABLE  2 


For  a  second  class  of  problems  we  consider  a  system  of  ordinary  differential 
equations 


x'  -  f(x,t),  xGRm  ,  teR1 


(15) 


which  is  assumed  to  have  a  2n  periodic  solution  Xg(t) .  If  x.  =  x^(t^) ,  i  e  N 
denotes  some  approximation  of  this  solution,  then  the  backward  Euler  discreti¬ 
zation 


xi+l  -  xi  +  hf<xl+l'ti+l) 


(16) 


represents  a  system  of  the  form  (5)  with  a  circle  mapping  (14a)  for  which 
w  -  1  .  Hence  we  can  apply  the  recursive  solution  process  to  this  system. 

As  a  first  sample  problem  we  consider  the  following  simple  form  of  van 
der  Pol's  equation 


x  -  A ( 1  -  x  )  x'  +  x  -  0 
As  shown,  for  example  in  [2],  the  solution  satisfies 


(16) 


x  -  2  cos(a>t)  +  A  (0.75  sin(ajt)  -  0.25  sin(3u>t)  +  0(A  ),  as  A  -*  0  (17) 

2 

and  oj-l  +  0(A)asA-»0. 
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In  Table  3  we  give  some  of  the  results  of  a  calculation  with  the  above 
mentioned  code  for  this  problem  when  w  -  1. 


n-64 

n-128 

n-256 

n-512 

T 

T 

T 

T 

1  ^ 

err 

nk 

err 

nk 

err 

nk 

err 

nk 

|  0.02 

0.01654 

3.91 

0.01560 

3.70 

0.01524 

3.61 

0.01511 

3.55  | 

|  0.04 

0.03358 

4.06 

0.03162 

3.71 

0.03081 

3.63 

0.03050 

3.54  | 

|  0.06 

0.05112 

4.06 

0.04811 

3.83 

0.04684 

3.63 

0.04631 

3.57  | 

|  0.08 

0.06917 

4.22 

0.06522 

3.75 

0.06345 

3.63 

0.06271 

3.57  | 

|  0.10 

0.08777 

4.22 

0.08292 

3.83 

0.08073 

3.67 

0.07979 

3.56  | 

TABLE  3 

Here  err  denotes  the  error  in  the  infinity  norm  of  Rn  between  the  com¬ 
puted  derivatives  x'  and  the  x' -values  obtained  from  the  expansion  (17) .  The 
corresponding  errors  in  the  x-values  are  slightly  smaller.  Evidently,  the 
error  due  to  the  zero-order  approximation  of  the  frequency  u>  dominates  the 
discretization  error.  Once  again  the  time  depends  linearly  on  n. 

As  a  second  problem  we  used  the  periodically  forced  brusselator 

2  2 
x'  -A+xy-Bx-x+a  coswt  ,  y’  -  Bx  -  x  y  (18) 

Here  a  and  w  are  parameters  and  a  frequent  choice  for  the  constants  A  and  B  is 

A  -  0.4,  B  —  1.2.  The  unforced  brusselator  (a  -  0)  then  has  a  limit  cycle  with 

the  natural  frequency  -  0.3750375  . 

For  our  computations  we  considered  the  resonance  case  u>  -  Wq.  Then  for 
0  <  | a |  <  a  ”  0.0033  three  periodic  solutions  are  obtained  two  of  which 
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are  inside  the  original  limit  cycle  and  one  outside.  With  increasing  |a|  the 
inside  solutions  approach  each  other  until  for  |a|  -  acrit  t*ley  coincide  and 
for  | a |  >  acrit  they  disappear.  The  amplitude  of  the  outside  solution 
increases  with  |a| 

In  order  to  obtain  suitable  starting  solutions  for  the  continuation  pro¬ 
cess  along  these  three  solution  branches,  we  imbed  (18)  into  a  homotopy 
depending  on  a  new  parameter  A.  For  this  the  equations  are  first  transformed 
to  the  new  time  t  ->  cj^t  .  Then,  by  adding  them,  we  see  that  z  =  x  +  y  satis¬ 
fies 

z'  -  T(A  -  x  -  a  cos  t)  ,  t  -  1A>q 

Now  x,x'  are  linear  expressions  of  z',z''  and  by  substituting  them  into  the 
transformed  equations  and  eliminating  x  and  y  ,  we  obtain  for 
u  -  z  -  (A  +  B/A)  the  equation 

u"  -  -f(u,u'  ,t,u>)  +  -u(AT-u')2  +  u'2(2A-A/B)  (19a) 

f(u,v,t,c<j)  -  T(1+A2-B)v  +  wv3  +  aT2cos  t  +  aT(AT-u' )  2sin  t]  (19b) 
From  this  the  desired  homotopy  is  constructed  as  follows 

u'  -  v,  v'  -  -Af (u,v, t,w(A) )  -  u(AT-A2v)2  +  A2v2 (2A-B/A) ,  (20) 
where  w(A)  -  A  -  0.0249625A.  Clearly,  for  A  -  1  this  system  is  equivalent  with 
(19a/b)  and  for  A  -  0  it  reduces  to  the  system  u’  -  v,  v'  -  -u  which  has  an 
infinite  set  of  2 it—  periodic  solutions.  The  limit  of  the  periodic  solutions  of 
(20)  for  A  -♦  0  can  be  derived  by  the  classical  small -parameter  technique.  For 
this  we  need  to  solve  the  equations 

J2,rf(u0,v0,  t,0)sin  t  dt  -  0,  J2wf (Uq  , Vq ,  t , 0)cos  ;  dt  -  0 


with 
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Uq  -  a  sin  t  +  p  cos  t,  v^  -  a  cos  t  -  0  sin  t 
For  a>(0)  -  A  the  evaluation  of  these  integrals  produces  the  system 

(1+A2— B)£  -  0.75A2£(q2  +  p2)  +  a  -  0,  (1+A2-B)a  -  0.75A2a(a2  +  p2)  +  a/A  =  0 

which  has  the  solutions 


a  -  Cy ,  0  - 


2  ,  B-l-A*  .  1/2 
-An,  C  -  -  ( - —) 

3(1+A  ) 


where  7  satisfies 


7  -  7  +  D  -  0 ,  D  - 


4a 


3(AC)3(1+A2) 


4/3 

This  equation  for  7  has  three  real  roots  for  |D|  <  2/3  ;  that  is,  for 


(21) 


|  a  |  <  a 

-  cnt 


2a 


34'3D 


»  0.0033 


Exactly  one  of  these  roots  has  a  different  sign  than  a  and  always  exceeds  1  in 

modulus.  The  other  two  roots  have  the  same  sign  as  a  and  are  less  than  one  in 

modulus;  moreover,  they  coincide  for  faf  -  acrit  and  become  complex  for 

| a|  >  a  .  . 

cnt 


For  the  computation  we  use  now  the  homotopy  (20)  with  a  -  0  in  which  case 
the  three  roots  of  (21)  are  7  -  -1,0,1  .  With  the  corresponding  functions 
Uq,  Vq  for  one  of  these  roots  as  starting  point  we  begin  the  continuation  pro¬ 
cess  from  A  -  0  and  proceed  until  at  A  -  1  some  functions  u^,v^  are 
obtained.  They,  in  turn,  become  the  starting  point  for  the  continuation  pro¬ 
cess  in  the  parameter  a  beginning  from  a  -  0  and  with  A  held  at  the  fixed 
value  1.  During  both  continuations  it  is  advisable  to  introduce  in  (21)  the 
transformed  variables  u  -»  u  -  u ^  ,  v  -*  v  -  v .  with  j-0,1,  respectively. 

In  terms  of  the  original  variables  x,y  of  (18)  the  second  continuation 
process  produces  points  along  a  branch  of  periodic  solutions  of  the  forced 
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brusselator.  The  properties  of  these  solutions  depend  on  the  value  of  7  with 
which  the  first  continuation  was  started.  More  specifically,  as  expected,  in 
the  two  cases  7  -  —1,0  these  solutions  cease  to  exist  when  |a|  exceeds  acr^t 
while  for  7  -  1  our  second  continuation  process  did  not  fail. 


We  present  here  only  some  results  for  the  case  latter  case  7  -  1.  The 
first  continuation  always  proceeded  with  steps  of  order  0.1  from  X  -  0  to 
X  -  1  .  In  Table  4  we  give  some  results  for  the  second  continuation  with 
n  —  64  . 


1  a 

r 

_L 

nk 

a 

r 

— 

nk 

|  0.00 

0.4561 

4.22 

0.11 

0.3938 

4.30  | 

|  0.01 

0.4292 

4.18 

0.12 

0.3865 

4.38  | 

|  0.02 

0.3981 

4.32 

0.13 

0.3579 

4.38  | 

|  0.03 

0.3604 

4.38 

0.14 

0.3363 

4.32  | 

|  0.04 

0.3119 

4.34 

0.15 

0.3216 

4.32  | 

|  0.05 

0.2399 

4.22 

0.20 

0.2960 

4.42  | 

)  0.06 

0.08461 

4.48 

0.25 

0.2971 

4.38  | 

j  0.07 

0.02140 

4.22 

0.30 

0.3070 

4.48  | 

|  0.08 

0.03430 

4.30 

0.35 

0.3224 

4.38  | 

|  0.09 

0.09779 

4.42 

0.40 

0.3424 

4.61  | 

|  0.10 

0.1986 

4.30 

TABLE  4 


2  2  1/2 

Here  r  -  (  x(0)  +  x'(0)  )  and  not  all  continuation  steps  are  shown. 


the  process  slows  down  before  the  turning  points  of  r 


In  fact, 


The  linear 


dependence  of  the  time  on  n  is  again  evident. 
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